NONPROFIT DATA

Create NPO Data Frame

print("Don't run me")

# Import subset from NCCS BMF Data on nonprofits in Syr MSA

source_data("https://github.com/lecy/analyzing-nonprofit-service-areas/blob/master/DATA/syr_orgs.Rda?raw=true")


# Subset - Extract poverty orgs

pov_orgs <- syr_orgs[which (  syr_orgs$NTEECC == "B61" | 
                              syr_orgs$NTEECC == "B82" |
                              syr_orgs$NTEECC == "B92" |
                            syr_orgs$nteeFinal1 == "E" | 
                            syr_orgs$nteeFinal1 == "F" |
                            syr_orgs$nteeFinal1 == "I" |
                              syr_orgs$NTEECC == "K30" |
                              syr_orgs$NTEECC == "K31" |
                              syr_orgs$NTEECC == "K34" |
                              syr_orgs$NTEECC == "K35" |
                              syr_orgs$NTEECC == "K36" |
                            syr_orgs$nteeFinal1 == "L" |  
                            syr_orgs$nteeFinal1 == "O" |
                            syr_orgs$nteeFinal1 == "P" |  
                            syr_orgs$nteeFinal1 == "R" |  
                            syr_orgs$nteeFinal1 == "S" |
                            syr_orgs$nteeFinal1 == "T"), ]

Geocode Poverty Organizations

print("Don't run me")

# download specific ggmap, register Google key

devtools::install_github("dkahle/ggmap")

register_google(key = 'AIzaSyAMrECdXqpOt443ms6nzh118uUsqC2lE6M',
                account_type = "premium", day_limit = 15000)

# subset orgs for gecoding address information

pov_orgs_gps <- subset(pov_orgs, select=c(EIN:INCOME))

pov_orgs_gps$whole_address <- do.call(paste, c(pov_orgs_gps[ 
                                      c("ADDRESS", "CITY","STATE",
                                        "ZIP5")], sep = ", ")) 

pov_orgs_gps <- subset(pov_orgs_gps, select=c(whole_address, EIN:INCOME))

for (i in 1:nrow(pov_orgs_gps)) {
  latlon = geocode(pov_orgs_gps[i,1]) 
  pov_orgs_gps$lon[i] = as.numeric(latlon[1])
  pov_orgs_gps$lat[i] = as.numeric(latlon[2])
}

# Make sure you find a way to then take any PO Boxes and place them in the 
# center of centroids, and only the ones that are NAs.

CENSUS DATA

Download & clean Census ACS 2011-2015 5-yr data

Census API Package

print("Don't run me")

# Census API 
# install.packages("devtools")
# devtools::install_github("hrecht/censusapi")

censuskey <- "f8064a17832b439e690be95ab0e0ef9a78617584"

var.list <- c("NAME", 
              "B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
              "B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
              "B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E",                      "B15003_017E", "B15003_018E",
              "B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
              "B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E", 
              "B03001_002E", "B03001_003E", 
              "GEOID" )

acs_2015 <- getCensus( name="acs5",
                       vintage = 2015,
                       key = censuskey,
                       vars= var.list, 
                       region ="tract:*", 
                       regionin ="state:36" 
                     )


# subset to Syracuse 
# sum(acs_2015$county == "067")
# sum(acs_2015$county == "075")
# sum(acs_2015$county == "053")



acs_2015_syr <- acs_2015[which (acs_2015$county == "067" | 
                                acs_2015$county == "075" |
                                acs_2015$county == "053"), ]

Rename Variables

print("Don't run me")

# Race variables are single race, non-hispanic. For example, "white" is really "white non-hispanic." This is true for every race category listed except, obviously, hispanic.


rename( acs_2015_syr, 
       c("B25077_001E" = "mdn_hous_val", 
         "B25003_001E" = "tenure_tot", 
         "B25003_002E" = "tenure_own", 
         "B25003_003E" = "tenure_rent",
         "B25002_001E" = "tot_occup", 
         "B25002_002E" = "occupied",
         "B25002_003E" = "vacant", 
         "B01003_001E" = "tot_pop", 
         "B19013_001E" = "mhh_income",
         "B17020_002E" = "poverty",
         "B12006_001E" = "labor_part",
         "B23025_005E" = "unemploy",
         "B15003_017E" = "high_sch",
         "B15003_018E" = "ged",
         "B03001_001E" = "tot_race",
         "B03002_003E" = "white",             
         "B03002_004E" = "black", 
         "B03002_005E" = "am_ind",
         "B03002_006E" = "asian", 
         "B03002_007E" = "islander",
         "B03002_008E" = "other", 
         "B03002_009E" = "mixed",
         "B03001_002E" = "not_hispanic",
         "B03001_003E" = "hispanic"
      
           ))

acs_2015_syr <- rename.vars(acs_2015_syr,
            c("B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
              "B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
              "B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E",                      "B15003_017E", "B15003_018E",
              "B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
              "B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E", 
              "B03001_002E", "B03001_003E"),  
            c("mdn_hous_val", "tenure_tot", "tenure_own", "tenure_rent",
              "tot_occup", "occupied", "vacant", "tot_pop", 
              "mhh_income", "poverty", "labor_part", "unemploy",
              "high_sch", "ged",
              "tot_race", "white", "black", "am_ind",
              "asian", "islander", "other", "mixed",
              "not_hispanic", "hispanic"))

# Create census race proportion variable

acs_2015_syr$porp_white <- (acs_2015_syr$white) / (acs_2015_syr$tot_race)
acs_2015_syr$porp_m <- abs((acs_2015_syr$porp_white) - 1) 

names( acs_2015_syr )

MAPS

Load in Census Shapefiles from GitHub

Note: To see how these shapefiles were created, check Shapefiles to Geojson Files HTML file for more information.

file.url <- "https://raw.githubusercontent.com/lecy/analyzing-nonprofit-service-areas/master/SHAPEFILES/syr_merged.geojson"

syr <- geojson_read( file.url, method="local", what="sp")

spTransform(syr, CRS("+proj=longlat +datum=WGS84"))
## class       : SpatialPolygonsDataFrame 
## features    : 186 
## extent      : -76.64471, -75.23954, 42.72384, 43.70702  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 68
## names       : STATEFP.1, COUNTYFP.1, TRACTCE.1,     GEOID.1, NAME.1,          NAMELSAD.1, MTFCC.1, FUNCSTAT.1,   ALAND.1, AWATER.1,  INTPTLAT.1,   INTPTLON.1, STATEFP.2, COUNTYFP.2, TRACTCE.2, ... 
## min values  :        36,        053,    030101, 36053030101, 301.01, Census Tract 301.01,   G5020,          S, 102677781,        0, +42.7960459, -075.3432945,        36,        067,    000100, ... 
## max values  :        36,        053,    031100, 36053031100,    311,    Census Tract 311,   G5020,          S,   8864612,    93115, +43.1261230, -075.8928470,        36,        067,    940000, ...
file.url <- "https://raw.githubusercontent.com/lecy/analyzing-nonprofit-service-areas/master/SHAPEFILES/syr_merged_cen.geojson"

centroids <- geojson_read( file.url, method="local", what="sp")

spTransform(centroids, CRS("+proj=longlat +datum=WGS84"))
## class       : SpatialPointsDataFrame 
## features    : 186 
## extent      : -76.55702, -75.34358, 42.78044, 43.63666  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : dummy 
## min values  :     0 
## max values  :     0
file.url <- "https://raw.githubusercontent.com/lecy/analyzing-nonprofit-service-areas/master/SHAPEFILES/pov_orgs.geojson"

pov_orgs <- geojson_read( file.url, method="local", what="sp")

spTransform(pov_orgs, CRS("+proj=longlat +datum=WGS84"))
## class       : SpatialPointsDataFrame 
## features    : 1095 
## extent      : -76.58932, -75.27822, 42.7584, 43.58445  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 21
## names       :                                 whole_address,       EIN,  FIPS, NTEECC, FILER, ZFILER,                               NAME,                  ADDRESS,          CITY, STATE,  ZIP5,  GEN, SUBSECCD, RULEDATE, FNDNCD, ... 
## min values  :      1 ARKIE ALBANESE AVE, MANLIUS, NY, 13104, 010587946, 36053,    B61,     N,      N, 0SWEGO COUNTY DEPUTIES ASSOCIATION,     1 ARKIE ALBANESE AVE, BALDWINSVILLE,    NY, 13027, 0000,       02,   000000,     00, ... 
## max values  : SUNY COLLEGE OF FORESTRY, SYRACUSE, NY, 13210, 952505709, 36075,    T99,     Y,      Y,       ZONTA INTERNATIONAL SYRACUSE, SUNY COLLEGE OF FORESTRY, WEST EDMESTON,    NY, 13485, 9434,       92,   201607,     21, ...
par( mar=c(0,0,0,0) )  # drop plot margins

Check shapefiles with plots

plot( syr )
plot( centroids, col="black", pch = 20, add = T)

plot( syr )
plot( centroids, col="black", pch = 20, add = T )
plot( pov_orgs, col="red", lwd = 2, add = T )

names(syr)
##  [1] "STATEFP.1"    "COUNTYFP.1"   "TRACTCE.1"    "GEOID.1"     
##  [5] "NAME.1"       "NAMELSAD.1"   "MTFCC.1"      "FUNCSTAT.1"  
##  [9] "ALAND.1"      "AWATER.1"     "INTPTLAT.1"   "INTPTLON.1"  
## [13] "STATEFP.2"    "COUNTYFP.2"   "TRACTCE.2"    "GEOID.2"     
## [17] "NAME.2"       "NAMELSAD.2"   "MTFCC.2"      "FUNCSTAT.2"  
## [21] "ALAND.2"      "AWATER.2"     "INTPTLAT.2"   "INTPTLON.2"  
## [25] "STATEFP"      "COUNTYFP"     "TRACTCE"      "GEOID"       
## [29] "NAME"         "NAMELSAD"     "MTFCC"        "FUNCSTAT"    
## [33] "ALAND"        "AWATER"       "INTPTLAT"     "INTPTLON"    
## [37] "GEOID_full"   "NAME.3"       "GEOID.3"      "state"       
## [41] "county"       "tract"        "mdn_hous_val" "tenure_tot"  
## [45] "tenure_own"   "tenure_rent"  "tot_occup"    "occupied"    
## [49] "vacant"       "tot_pop"      "mhh_income"   "poverty"     
## [53] "labor_part"   "unemploy"     "high_sch"     "ged"         
## [57] "tot_race"     "white"        "black"        "am_ind"      
## [61] "asian"        "islander"     "other"        "mixed"       
## [65] "not_hispanic" "hispanic"     "porp_white"   "porp_m"

Create zoom enabled map of Syracuse, NY

map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY")

map %>% addProviderTiles(providers$CartoDB)

Set the view of the Syracuse MSA

syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr, fillOpacity = 0.3) %>%
  addCircles(data = pov_orgs) %>%
  addCircles(data = centroids)

syr_map %>% addProviderTiles(providers$CartoDB)  

Edit color format of polygons and data points

syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr, fillOpacity = 0.3, weight = 2) %>%
  addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
  addCircles(data = centroids, color= "black", opacity = 1.0)

syr_map %>% addProviderTiles(providers$CartoDB)  

Adding census tract information to map

# change variables to be numeric

syr$porp_white <- as.numeric(as.character(syr$porp_white))
syr$porp_m <- as.numeric(as.character(syr$porp_m))
syr$mhh_income <- as.numeric(as.character(syr$mhh_income))
class(syr$porp_m)
## [1] "numeric"
class(syr$porp_white)
## [1] "numeric"
class(syr$mhh_income)
## [1] "numeric"
# By race

pal <- colorNumeric(
  palette = "Blues",
  domain = syr$porp_m )


syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr, fillOpacity = 0.3, weight = 2,
              color = ~pal(syr$porp_m)) %>%
  addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
  addCircles(data = centroids, color= "black", opacity = 1.0) 


syr_map %>% addProviderTiles(providers$CartoDB)
# By Income level

pal2 <- colorNumeric(
  palette = "Blues",
  domain = syr$mhh_income )


syr_map2 <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr, fillOpacity = 0.3, weight = 2,
              color = ~pal2(syr$mhh_income)) %>%
  addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
  addCircles(data = centroids, color= "black", opacity = 1.0) 


syr_map2 %>% addProviderTiles(providers$CartoDB)

ANALYSIS

Create Service Catchment Areas

print("Don't run me")


# Create buffer, calculate number of nonprofit points around centroids

centroids_buff <- gBuffer(spgeom = centroids, width = .024)

spTransform(centroids_buff, CRS("+proj=longlat +datum=WGS84"))

CRS.new<-CRS("+proj=longlat +datum=WGS84")
proj4string(centroids_buff) <- CRS.new 
proj4string(pov_orgs) <- CRS.new

over <- over( pov_orgs, centroids_buff)
# syr$locate <- over$GEOID_full
# syr$count <- rep.int(1, 186)

# other option?
# syr$test <- 1:nrow(over_clean)



pal2 <- colorNumeric(
  palette = "Blues",
  domain = syr$mhh_income )

syr_map2 <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr, fillOpacity = 0.3, weight = 2,
              color = ~pal2(syr$mhh_income)) %>%
  addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
  addCircles(data = centroids, color= "black", opacity = 1.0) %>%
  addPolygons(data = centroids_buff, fillOpacity = 0.3, weight = 2,)
             
syr_map2 %>% addProviderTiles(providers$CartoDB)

# Not quite sure what is needed. How do I add up values for every centroid from the overlay?

# Similar to above, but combines over and removing NAs into one step
# pov_orgs_t <- pov_orgs 
# pov_orgs_t$test <- pov_orgs[centroids_buff,]

Analysis continued

# Poportion of census tract that are high/low diversity
  # number of nonprofits by type in high/low racially diverse areas
  # number of nonprofits by type in low/high income census tracts


# Corr analysis between nonprofit access and census tract attributes
  # Number of nonprofits location by census tract income or race
    # Race
    # Income

Final section: Account for all urban census tracts in 25 cities